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Abstract: The problem of effectively combining data with a mathematical model constitutes a major 
challenge in applied mathematics. It is particular challenging for high-dimensional dynamical systems 
where data is received sequentially in time and the objective is to estimate the system state in an 
on-line fashion; this situation arises, for example, in weather forecasting. The sequential particle filter 
is then impractical and ad hoc filters, which employ some form of Gaussian approximation, are widely 
used. Prototypical of these ad hoc filters is the 3DVAR method. The goal of this paper is to analyze 
the 3DVAR method, using the Lorenz '63 model to exemplify the key ideas. The situation where 
the data is partial and noisy is studied, and both discrete time and continuous time data streams 
are considered. The theory demonstrates how the widely used technique of variance infiation acts to 
stabilize the filter, and hence leads to asymptotic accuracy. 

1. Introduction 

Data assimilation concerns estimation of the state of a dynamical system by combining observed data with 
the underlying mathematical model. It finds widespread application in the geophysical sciences, including 
meteorology [16], oceanography [2] and oil reservoir simulation [23]. Both filtering methods, which update 
the state sequentially, and variational methods, which can use an entire time window of data, are used [1]. 
However, the dimensions of the systems arising in the applications of interest are enormous - of O(IO^) in 
global weather forecasting, for example. This makes rigorous Bayesian approaches such as the sequential 
particle filter [6] , for the filtering problem, or MCMC methods for the variational problem [28] , prohibitively 
expensive in on-line scenarios. 

For this reason various ad hoc methodologies are typically used. In the context of filtering these usually rely 
on making some form of Gaussian ansatz [31]. The 3DVAR method [17, 25] is the simplest Gaussian filter, 
relying on fixed (with respect to the data time-index increment) forecast and analysis model covariances, 
related through a Kalman update. A more sophisticated idea is to update the forecast covariance via the 
linearized dynamics, again computing the analysis covariance via a Kalman update, leading to the extended 
Kalman filter [14]. In high dimensions computing the full linearized dynamics is not practical. For this 
reason the ensemble Kalman filter [7, 8] is widely used, in which the forecast covariance is estimated from 
an ensemble of particles, and each particle is updated in Kalman fashion. An active current area of research 
in filtering concerns the development of methods which retain the computational expediency of approximate 
Gaussian filters, but which incorporate physically motivated structure into the forecast and analysis steps 
[21, 20], and are non-Gaussian. 

Despite the widespread use of these many variants on approximate Gaussian filters, systematic mathemat- 
ical analysis remains in its infancy. Because the 3DVAR method is prototypical of other more sophisticated 
ad hoc filters it is natural to develop a thorough understanding of the mathematical properties of this filter. 
Two recent papers address these issues in the context of the Navier-Stokes equation, for data streams which 
are discrete in time [5] and continuous in time [4]. These papers study the situation where the observations 
are partial (only low frequency spatial information is observed) and subject to small noise. Conditions are 
established under which the filter can recover from an order one initial error and, after enough time has 
elapsed, estimate the entire system state to within an accuracy level determined by the observational noise 
scale; this is termed filter accuracy. Key to understanding, and proving, these results on the 3DVAR filter 
for the Navier-Stokes equation are a pair of papers by Titi and co-workers which study the synchronization 
of the Navier-Stokes equation with a true signal which is fed into only the low frequency spatial modes of the 
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system, without noise [24, 13]; the higher modes then synchronize because of the underlying dynamics. The 
idea that a finite ammount of information effectively governs the large-time behaviour of the Navier-Stokes 
equation goes back to early studies of the equation as a dynamical system [10] and is known as the determin- 
ing node or mode property in the modern literature [26]. The papers [5, 4] demonstrate that the technique 
of variance inflation, widely employed by practitioners in high dimensional filtering, can be understood as a 
method to add greater weight to the data, thereby allowing the synchronization effect to take hold. 

The Lorenz '63 model [18, 27] provides a useful metaphor for various aspects of the Navier-Stokes equation, 
being dissipative with a quadratic energy-conserving nonlinearity [9]. In particular, the Lorenz model exhibits 
a form of synchronization analogous to that mentioned above for the Navier-Stokes equation [13]. This 
strongly suggests that results proved for 3DVAR applied to the Navier-Stokes equation will have analogies 
for the Lorenz equations. The purpose of this paper is to substantiate this assertion. 

The presentation is organized as follows. In section 2 we describe the Bayesian formulation of the inverse 
problem of sequential data assimilation; we also present a brief introduction to the relevant properties of 
the Lorenz '63 model and describe the 3DVAR filtering schemes for both discrete and continuous time data 
streams. In section 3 we derive Theorem 3.2 concerning the 3DVAR algorithm applied to the Lorenz model 
with discrete time data. This is analogous to Theorem 3.3 in [5] for the Navier-Stokes equation. However, in 
contrast to that paper, we study Gaussian (and hence unbounded) observational noise and, as a consequence, 
our results are proved in mean square rather than almost surely. In section 4 we extend the accuracy result 
to the continuous time data stream setting: Theorem 4.1; the result is analogous to Theorem 4.3 in [4] which 
concerns the Navier-Stokes equation. Section 5 contains numerical results which illustrate the theory. We 
make concluding remarks in section 6. 

2. Set-Up 

In subsection 2.1 we formulate the probabilistic inverse problem which arises from attempting to estimate 
the state of a dynamical system subject to uncertain initial condition, and given partial, noisy observations. 
Subsection 2.2 introduces the Lorenz '63 model which we employ throughout this paper. In subsections 2.3 
and 2.4 we describe the discrete and continuous 3DVAR filters whose properties we study in subsequent 
sections. 



2.1. Inverse Problem 

Consider a model whose dynamics is governed by the equation 

with initial condition u(0) ~ uq £ W. We assume the the initial condition is uncertain and only its statistical 
distribution is known, namely the Gaussian uq iV(mo, Cq). Assuming that the equation has a solution for 
any uq G MP and all positive times, we let ^'(•, •) : x M+ — )• M.P be the solution operator for equation (2.1). 
Now suppose that we observe the system at equally spaced times tk = kh for all k G Z+. For simplicity we 
write ^'(•) :— ^'(•; h). Defining Uk = u{tk) — ^(uq; kh) we have 

"fe+i = *("fc), k e Z+. (2.2) 

We assume that the data {yk}kei,+ is found from noisily observing a linear operator H applied to the system 
state, at each time tk, so that 

yk+i = Huk+i + Vk+i, fc e N. (2.3) 

Here {i'k}k&i is an i.i.d. sequence of random variables, independent of uq, with vi ^ N{0, T) and H denotes 
a linear operator from Rp to K™, with m < p. If the rank of H is less than p the system is said to be partially 
observed. We denote the accumulated data up to time k by Yk := {yj}^^i- The pair {uk,Yk) is a jointly 
varying random variable in W x K'^'™ . The goal of filtering is to determine the distribution of the conditioned 
random variable Uk\Yk, and to update it sequentially as k is incremented. This is an inverse problem for the 
system state, given observed data, and it has been regularized by means of the Bayesian formulation. 
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2.2. Forward Model: Lorenz '63 



When analyzing the 3DVAR approach to the filtering problem we will focus our attention on a particular 
model problem, namely the classical Lorenz '63 system [18]. In this section we introduce the model and 
summarize the properties relevant to this paper. The Lorenz equations are a system of three coupled non- 
linear ordinary differential equations whose solution w G M^, where u = {ux,Uy,Uz), satisfies 



a{uy-U3;), (2.4a) 
My - ~aux - Uy - u^Uz, (2.4b) 



Uy 

Uz — u^Uy ~ bUz ~ b{r + a). (2.4c) 
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Note that we have employed a coordinate system where origin is shifted to the point (0,0, — (r + a)) as 
discussed in [29]. Throughout this paper we will use the classical parameter values {a,b,r) = (10, |,28) in 
all of our numerical experiments. At these values, the system is chaotic [30] and has one positive and one 
negative Lyapunov exponent and the third is zero, reflecting time translation-invariance. Our theoretical 
results, however, simply require that a,b > 1 and r > and we make this assumption, without further 
comment, throughout the remainder of the paper. 

In the following it is helpful to write the Lorenz equation in the following form as given in [9], [12]: 

— +Au + B{u,u)^f, u{0)=uo, (2.5) 

where 

A = 

\ b J 

( ' 

B{u,u)=\ {u^Uz + UzUx)/2 

\ -{UxUy + UyUx)l2 

We use the notation (•, •) and | • | for the standard Euclidean inner-product and norm. When describing our 
observations it will also be useful to employ the projection matrices P and Q defined by 

(2.6) 

\ 1 y 

We will use the following properties of A and B: 
Properties 2.1 ([12]). 

1. {Au, u) — au^ + Uy + buj. > \u\^ provided that a,b > 1. 

2. {B{u,u),u) = 0. 

3. B{u,u) = B{U,u). 

I \B{u,u)\<2~^\u\\u\. 

5. |(B(u,u),u)| < 2~^\u\\u\\Pil\. 

We will also use the following: 

Proposition 2.2. ([12], Theorem 2.2) Equation (2.5) has a global attractor A. Let u be a trajectory with 
uq e A. Then \u(t)\'^ < K for allteR where 

Figure 1 illustrates the properties of the equation. Sub-figure 1(a) shows the global attractor A. Sub-figures 
1(b), 1(c) and 1(d) show the components u^, Uy and Uz^ respectively, plotted against time. 

3 



















• 





• 


Q = 


• 


1 


• 



















(a) Sample trajectory for a solution lying on global attractor. (b) component 




(c) Uy component. (d) Uz component. 

Fig 1. Lorenz attractor and individual components. 



2.3. 3DVAR: Discrete Time Data 

In this section we describe the 3DVAR filtering scheme for the model (2.1) in the case where the system is 
observed discretely at equally spaced time points. The system state at time ~ kh is denoted by Uk — u{tk) 
and the data upto that time is Yf. = {yj}j=i- Recall that our aim is to find the probability distribution of 
UfejYfc. Approximate Gaussian filters, of which 3DVAR is a prototype, impose the following approximation: 

F{uk\Yk)^N{mk,Ck). (2.8) 

Given this assumption the filtering scheme can be written as an update rule 

(mfc, Cfc) i-^ (m^+i, Cfe+i). (2.9) 

To determine this update we make a further Gaussian approximation, namely that Uk+i given follows a 
Gaussian distribution: 

V{uk+i\Yk) = N{mk+i,Ck+i). (2.10) 

Now we can break the update rule into two steps of prediction {mk,Ck) ^ (m^+i, C/t+i) and analysis 
(TOfe+i,Cfe+i) — > {mk+i,Ck+i)- For the prediction step we assume that m/c+i = '^{mk) whilst the choice 
of the covariance matrix Ck+i depends upon the choice of particular approximate Gaussian filter under 
consideration. For the analysis step, (2.10) together with the fact that yk+i\uk+i ^ N{Huk+i,T) and 
application of Bayes' rule, implies that 

Ufe+i|Yfc+i ~ N{mk+i,Ck+i) (2-11) 

where [11] 

Ck+i = Ck+i — Ck+iH* [T + HCk+iH*y^HCk+i (2.12a) 
ruk+i = 'f{mk) + Ck+iH*{r + HCk+iH*)-\yk+i-H^{mk)). (2.12b) 



As mentioned the choice of update rule Ck Ck+i defines the particular approximate Gaussian filtering 
scheme. For the 3DVAR scheme we impose Ck+i = C Vfc G N where C is a positive definite p x p matrix. 
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From equation (2.12b) we then get 

ruk+i = ^imk) + CH*{T + HCH*)-\y„+,-H'^{mk)) 

= {I -GH)^{mk)+Gyk+i (2.13) 

where 

G:^CH*{T + HCH*y^ (2.14) 

is called Kalman gain matrix. Equation (2.13) is analyzed in section 3. 

Another way of defining the 3DVAR filter is by means of the following variational definition: 

/ \ 1 2 1 1 2\ 

mfc+i=argmin (-\\C~i{m-^{mk))\\ + -\\T-^yk+i - Hm)\\ j. (2.15) 

This coincides with the previous definition because the mean of a Gaussian can be characterized as the 
minimizer of the negative of the logarithm of the probability density function. That this negative logarithm 
is the sum of two squares follows from Bayes' theorem. From the variational formulation, it is clear that the 
3DVAR filter is a compromise between fitting the model and the data. The model uncertainty is characterized 
by a fixed covariance C, and the data uncertainty by a fixed covariance F; the ratio of the size of these two 
covariances will play an important role in what follows. 



2.4. 3DVAR: Continuous Time Data 

In this section we describe the limit of high frequency observations h which, with appropriate scaling 
of the noise covariance with respect to the observation time h, leads to a stochastic differential equation 
(SDK) Umit for the 3DVAR fiher. We refer to this SDE as the continuous time 3DVAR filter. We give a brief 
derivation, referring to [4] for further details and to [3] for a related analysis of continuous time limits in the 
context of the ensemble Kalman filter. 

We assume the following scaling for the observation error covariance matrix: F = ji^o- Thus, although the 
data arrives more and more frequently, as we consider the limit /i — 0, it is also becoming more uncertain; this 
trade-off leads to the SDE limit. Define the sequence of variables {zk}keN by the relation Zk+i = zj- + hy^+i 
and Zq — 0. Then 

Zk+i ^ Zk + hHuk+i + \/hrojk, zo = 0. (2-16) 
Here jk N{0, 1). By rearranging and taking limit as /i — >■ we get 

where w is an K™ valued standard Brownian motion. We think of Z{t) := {z(s)}sg[o,t] as being the data. For 
each fixed t we have the jointly varying random variable (u(i), Z{t)) gRp x C([0, t];MJ"). We are interested 
in the filtering problem of determining the sequence of conditioned probability distributions implied by the 
random variable u{t)\Z{t) in M.^. The 3DVAR filter imposes Gaussian approximations of the form N{m{t), C) . 
We now derive the evolution equation for m(t). 

Recall the vector field J" which drives equation (2.1). Using equation (2.16) in (2.13), together with the 
fact that ^{u) = u + hF{u) + Oih"^), gives 

=mn + hT{mr,)+0{h^) + hCH*{ro + hHCH*)-^ ^ z^+i^ z,, _^ jj^^ ^ (2.I8) 

Rearranging and taking limit /i — !■ gives 

^ = Hm) + CH*r^' (^g - iJm) . (2.19) 

Equation (2.19) defines the continuous time 3DVAR filtering scheme and is analyzed in section 4. The data 
should be viewed as the continuous time stream Z{t) — {z(s)}sg[o,t] and equations (2.17) and (2.19) as 
stochastic differential equations driven by w and z respectively. 
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3. Analysis of Discrete Time 3DVAR 



In this section we analyse the discrete time 3DVAR algorithm when applied to a partially observed Lorenz 
'63 model; in particular we assume only that the Uj; component is observed. We start, in subsection 3.1, 
with some general discussion of error propagation properties of the filter. In subsection 3.2 we study mean 
square behaviour of the filter for Gaussian noise. Recall the projection matrices P and Q given by (2.6), we 
will use these in the following. We will also use {vk} to denote the exact solution sequence from the Lorenz 
equations which underlies the data; this is to be contrasted with {uk} which denotes the random variable 
which, when conditioned on the data, is approximated by the 3DVAR filter. 



3.1. Preliminary Calculations 



Throughout we assume that the model covariance C = a^I and that H = (1, 0, 0), so that only Ux is observed 
may be written 



We also assume that T ~ e . Then the Kalman gain matrix is G = ^^H* , where n := The 3DVAR filter 



ruk+i = [j^P + Q] *(™fe) + -^Vk+iH*. (3.1) 

We define v to be the true solution of the Lorenz equation (2.5) which underlies the data, and we define 
Vk = v{kh), the solution at observation times. Note that, since F = e^, it is consistent to assume that the 
observation errors have the form 

/ 

I'fc = 
V 

where are i.i.d. random variables on M. We will consider the case ~ -/V(0, 1) for simplicity of exposition. 
Note that we may write 

Vk+iH* = Pvk+i + Vk+i 
= P^{vk) + vu+i. 



Thus 



Observe that 



mk+i = ( -^P + q] *(mfe) + (P^{vk) + lyk+i) ■ (3.2) 



Vk+i= (j^P + Q)-^{vk) + j^P'^{vk)- (3.3) 

We are interested in comparing nik, the output of the filter, with Vk the true signal which underlies the 
data. We define the error process S{t) as follows: 



m 



^ l^it) '^t t = tk 

■^{mk,t-tk)-v{t) if t e (tfc,tfc+i) 



Observe that 5 is discontinuous at times tj which are multiples of h, since rrifc+i 7^ ^[^k] h). In the following 
we write 5{tJ) for Ynn^^^- 6{t) and we define 5j — 5{tj). Thus 6j 7^ ^{tj)- Subtracting (3.3) from (3.2) we 
obtain 

Now consider the time interval {tk,tk+i)- Since 5{t) is simply given by the difference of two solutions of the 
Lorenz equations in this interval, we have 

— + AS + B{v,S) + B{S,v) + B{6,S)^0, t e {tk,tk+i). (3.5) 
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Taking the Euclidean inner product of equation (3.5) with S gives 



+{AS,d) + {B{v,S),S) + {B{5,v),6) + {B{6,S),6)^0 (3.6) 



2 dt 

which, on simphfying and using Properties 2.1, gives 



+{Ad,6)+2{B{v,S),6)^0, (3.7) 



2 dt 

and hence ^ 

l^ + \Sf + '2{B{v,S),5)<0. (3.8) 

In order to use (3.4) we wish to estimate the behaviour of <5(ijJ^]^) in terms of 6k- The foUowing is useful 
in this regard and may be proved by using (3.8) together with Properties 2.1(4). Note that K is defined by 
equation (2.7) and is necessarily greater than or equal to one, since 6, a > 1. 

Proposition 3.1 ([12]). Assume the true solution v lies on the global attractor A so that sup(>g|u(t)p < K 
with 

4(6-1) • 

Then for 13 = 2 (if^/^ i) it follows that \5{t)\^ < ISkfe^^^-^"^ forte [tk,tk+i). 



3.2. Accuracy Theorem 



In this subsection we assume that ^ ^(0, 1) and we study the behaviour of the filter in forward time 
when the size of the observational noise, ©(e), is small. The following result shows that, provided variance 
inflation is employed (77 small enough), the 3DVAR filter can recover from an 0(1) initial error and enter 
an 0{e) neighbourhood of the true signal. The results are proved in mean square. The reader will observe 
that the bound on the error behaves poorly as the observation time h goes to zero, a result of the over- 
weighting of observed data which is fluctuating wildly as /i — 0. This effect is removed in section 4 where 
the observational noise is scaled appropriately, in terms of ft. — >■ 0, to avoid this effect. 

For this theorem we define a norm || • || by = |up + |-Pup, where | • | is the Euclidean norm. 

Theorem 3.2. Let v be a solution of the Lorenz equation (2.5) with i;(0) G A, the global attractor. Assume 
that ^1 N{0, 1) so that the observational noise is Gaussian. Then there exist he > 0, A > such that for 
all rj sufficiently small and all h e (0, he) 

E||<5fc+i||' < (l-Aft)E||4||' + 2e2. (3.9) 

Consequently 

9,2 

limsupE||(5fe||^ < — . (3.10) 
fe->-oo Ara 

Proof. Recall that we have Ei^fc+i = and E|z^fe+i|^ = e^. On application of the projection P to the error 
equation (3.4) for 3DVAR we obtain 

E|P4+i|^ < (TT^)'E|i^^(t.-+i)l'+ (j^)'^^. (3.11) 

Since E|Q4+i|2 = ¥.\Q5{tl^^)\^ < E\6{t'^:^^)f we also obtain the bound 

m+i\' < l^^yE\PS{t-^,)f+mt-^,)f+(^^ye^ (3.12) 
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Define Mi and M2 by 



and 



M2(r) = -^(e-"--e-)+2f-^) e'"^ (3.14) 



Adding (3.11) to (3.12) and using Lemma 3.3 shows that 



1 



so that 



E|14+if <Mi(/i)E|4r + Af2(/i)E|P4|' + 2(^— j (3.15) 

E| 14+1 1 1' < A^(/i)E||4||' + ^j^, (3.16) 
where 

M(t) = max{Afi(T),M2(T)}. (3.17) 

Now we observe that 

A/i(0) = 1, A/i'(0) = -l + 2a( and 1/3(0) = 2^ — 

Vl + r?/ + 

Thus there exists a.n he > and a A > such that, for aU rj sufficiently small 

M(t,7?) < 1- At, Vre(0,/z,]. 

Hence the theorem is proved. □ 

The following lemma is used in the preceding proof. 
Lemma 3.3. Under the conditions of Theorem 3.2 for t G [tk,tk+i) we have 

|P5(t)|2 < ( e/3(*-*.) _ e-o(*-t.)) + |F4|2e-"(*-*'=) (3.I8) 



/3 + a 

and 



(3.19) 



,2 ^/VOIOfel 

/3 + a V + 1 1-a 

1 — a V / 
Proof. Taking inner product of (3.5) with P5, instead of with S as previously, we get 

1^ + {AS,PS)^0. (3.20) 

Let S = {Sx,Sy,Sz)'^ . Notice that \PSf — |4|^ and {AS, PS) = aS"^ — aS^Sy. Therefore equation (3.20) 
becomes 



By rearranging and applying Proposition 3.1 we get 



2 

- + a\P5\' < a\S{tk)\"e^^'-*'>. (3.21) 



dt 



Multiplying by integrating factor e"^* and integrating from t^. to t gives equation (3.18). 
Analysing the non-linear term in equation (3.8) with Property 2.1(5) gives 

\2{Biv,6),6)\ < \v\\PS\\S\ (3.22) 

< k'2\PS\\6\ (3.23) 

< ^K\PSf + ^\5f. (3.24) 

Substituting (3.18) and (3.24) in (3.8) gives 

^ + < ^^^M. fe^(*-*.) _ e-(*-t.)) + X|P4|2e-"(*-*^). (3.25) 
dt p + a \ ) 

Multiplying by the integrating factor e*^*"**"' and integrating from to t gives 

|J(t)|V-*^)~|4f < [ —^, - 1 U^^E^ U^-o.)it-t.) _ 1 . (3.26) 

p + a \ P + 1 1 — a / 1 — a V / 

Rearranging the above equation gives (3.19). □ 
4. Analysis of Continuous Time 3DVAR 

In this section we analyse application of the 3DVAR continuous filtering algorithm for the Lorenz equation 
(2.5). We will use (^)}te[o.oo) to denote the exact solution sequence from the Lorenz equations which 
underlies the data; this is to be contrasted with {u(t)}(g[o,oo) which denotes the random variable which, 
when conditioned on the data, is approximated by the 3DVAR filter. 

We study the continuous time 3DVAR filter in the case where = (1, 0, 0), Fo = e^, C = a^l and rj — 
To analyse the filter it is useful to have the truth v which gives rise to the data appearing in the filter itself. 
Thus (2.17) gives 

dz I — Aw 

We then eliminate z in equation (2.19) by using (4.1) to obtain 

— = T{m) + CH*Tg'H{v - m) + CH*T^ ' — . (4.2) 

In the specific case of the Lorenz equation we get 

— = -Am-B{m,m) + f + CH*T-^H{v-m) + CH*T^~^—. (4.3) 

From equation (4.2) with the choices of C, H and Fq detailed above we get 

d?Tz 1 6 Auj 

— — = -Am - B{m, m) + / + -P{v - to) + -P— (4.4) 
dt T] rj dt 

where we have extended w from a scalar Brownian motion to an M'^-valued Brownian motion for notational 
convenience. This has a unique global strong solution to G C([0, oo); M^). Indeed similar techniques used to 
prove the following result may be used to establish this global existence result, by applying the Ito formula 
to |mp and using the global existence theory in [22]; we omit the details. Recall K given by (2.7). 

9 



Theorem 4.1. Let m solve equation(4-4) o'f^ ^ solve equation (2.5) with initial data f (0) G A, the global 
attractor, so that sup^Q < K. Then for rjK < 4 we obtain 



E\m{t) - v{t)\' < e-^*|m(0) - ^^(0)|' + 4^(1 " e"^*)' (4-5) 



A = 2(1-^L (4.6) 



where X is defined by 



Thus 

2 

limsup4^^E|m(t) - v{t)\ < j-^. 
Proof. The true solution follows the model 

^ -Av - B(v,v) + f +-P(v ~v), (4.7) 
at ' rj 

where we include the last term, which is identically zero, for clear comparison with the filter equation (4.4). 
Define 5 = m — v and subtract equation (4.7) from equation(4.4) to obtain 

d(5 dz/^ 

— = -Am- B{m,m) + Av + B{v,v)-ri^^P5 + eri''^P — (4.8) 
at at 

= -AS-2B{v,S)- B(S,S)-r]-^PS + err^P—. (4.9) 

at 

Using Ito's formula gives 

ld\S\^ + {AS + 2B{v, S) + B{S, 6) + ^PS, S)dt < {cTj-^Pdw, S) + -^Tr (e^v^^P) dt. (4.10) 
2 rj 2 

Using Lemma 4.2 and the definition of A gives 

^d\5\^ + ^ \5\^dt < {er)-^Pdw, S) + ^Tt {e^v'^P) dt. (4.11) 
Rearranging and taking expectations gives 

(4.12) 

Use of the Gronwall inequality gives the desired result. □ 

The following lemma is used in the preceding proof. 
Lemma 4.2. Let v E A. Then 

{A6 + 2B{v, 5) + B{5, 5) + ^P5, 5) > (\ - \5\^ . (4.13) 
Proof. On expanding the inner product 

(AS + 2Biv, 5) + B{S, S) + -PS, 5) = (AS, S) + 2{B(v. S), S) + {B(S, 5),S) + (rj-'^PS, 5). (4.14) 

V 

We now use the Properties 2. 1(1), (5) and the fact that true solution lies on the global attractor so that 
\v\ < K. As a consequence we obtain 

(AS + 2B{v,S)+B(5,5) + -PS,S) > \Sf - KhsllPSl + -iPSf. (4.15) 
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Using Young's inequality with parameter 9 

{A6 + 2B{v,S) + B{6,S) + -PS,S) > |<5|' - ^^l^'^l' " + "l-P'^l'- (4-16) 

7y 29 Z rj 

Taking 9 — ^ yields the desired result 

{A6 + 2B{v, 5) + B{S, 6) + ^P5, S) > (^1 - \S\\ (4.17) 

□ 

5. Numerical Results 

In this section we present numerical results illustrating Theorems 3.2 and 4.1 established in the two preceding 
sections. All experiments are conducted with the parameters {a,b,r) = (10, |,28). Both the theorems are 
mean square results. However, our numerics are all based on a single realization of the filters in question. 

5.1. Discrete case 

Under the assumptions of Theorem 3.2 we expect the mean square error in \v — m\ to decrease exponentially 
until it is of the size of the observational noise squared. Hence we expect the estimate m to converge to a 
neighbourhood of the true solution v, where the size of the neighbourhood scales as the size of the noise 
which pollutes in observation, in mean square. The following experiment indicates that similar behaviour is 
in fact observed path wise. We set up the numerical experiments by computing the true solution v of the 
Lorenz equations using a fourth order Runge-Kutta method, and then adding Gaussian random noise to the 
observed x-component to create the data. We fix the parameter = 0.1 and vary the size of the observational 
noise e. 

We see in Figure 2 that initially the error |u(0) — to(0)| is of 0(1). It then decays exponentially with time 
and converges to 0(e). For this particular case we chose e — 0.316. A consequence of the second part of 
Theorem 3.2 is that the logarithm of the asymptotic mean squared error logE||(S|p varies linearly with the 
logarithm of the standard deviation of noise in the observations (e) . To compute the asymptotic mean square 
error we time-average the mean square error from a single long trajectory, until it reaches equilibrium. 
In Figure 3 we observe the log-linear decrease in the asymptotic error as the size of the noise decreases. 
Furthermore, the slope of the graph is approximately 2 as predicted by (3.10). 

5.2. Continuous case 

In the case of continuous observations we again compute a true trajectory of the Lorenz equation using 
a fourth order Runge-Kutta scheme. We then simulate the SDE (4.3) using the Euler-Maruyama scheme.^ 
Similarly to the discrete case, we consider path wise illustration of the mean square results in Theorem 4.1. In 
this case we fix the parameter rj = 0.1 and vary the size of observational error e. The initial error is expected 
to decay exponentially towards something of order e, and this is illustrated in Figure 4. Furthermore, the 
asymptotic mean square error is expected to be of order e^. As in the discrete case, we estimate this through 
time averaging and, in Figure 5, we observe the log-linear decrease in the asymptotic error as the size of the 
noise decreases. The slope of the graph is approximately 2, as predicted by equation (4.5). 



^Note that this is equivalent to creating the data z from (4.1) and solving (4.2) and, since we have access to the truth, is 
computationally expedient. 
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Fig 2. Decay of initial error from 0(1) to 0{e). 



Logarithm of asymptotic MSE 



logE||5|| 



2 
1 

■1 

2 
■2 

3 

■4 

-2.5 



-1.5 



-0.5 
log(e) 



0.5 



1.5 



Fig 3. Linear dependence of asymptotic logE||(5|p on loge for discrete observations. 
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Fig 4. Decay of initial error from 0(1) to 0{e). 
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Fig 5. Linear dependence of asymptotic logE|<5p on loge for discrete observations. 
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6. Conclusion 

The study of approximate Gaussian filters for the incorporation of data into high dimensional dynamical sys- 
tems provides a rich field for applied mathematicians. Potentially such analysis can shed light on algorithms 
currently in use, whilst also suggesting methods for the improvement of those algorithms. However, rigor- 
ous analysis of these filters is in its infancy. The current work demonstrates the properties of the 3DVAR 
algorithm when applied to the partially observed Lorenz '63 model; it is analogous to the more involved 
theory developed for the 3DVAR filter applied to the partially observed Navier-Stokes equations in [5, 4]. 
Work of this type can be built upon in three primary directions: firstly to consider other model dynamical 
systems of interest to practitioners, such as the Lorenz '96 model [19]; secondly to consider other observation 
models, such as pointwise velocity field measurements or Lagrangian data for the Navier-Stokes equations, 
building on the theory of determining modes [15]; and finally to consider more sophisticated filters such as 
the extended [14] and ensemble [7, 8] Kalman filters. 
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